Line Vortex Simulations

APAM1601y
Columbia University

Introduction

References: 
“Vortex crystals from 2D Euler flow; Experiment and simulation,” Phys. Fluids, 11(4), p . 905 (1999). 

See wonderful research at the UCSD Nonneutral Plasma Group.

This notebook computes the motion of a (very small) collection of line vortices using a "particle-in-cell" (PIC) simulation. This technique is common and found frequently in simulations of ionized gas (i.e. plasma). The key to a PIC simulation is the solution of the potentials on a fixed grid and the interpolation on this grid to find the local motion of the line vortices.

Equations

The equations of motion for the ith vortex line are
	(∂_t) (Subscript[x, i]) =  –(∂_y)Ψ((Subscript[x, i]), (Subscript[y, i]), t)
	(∂_t) (Subscript[y, i]) =   (∂_x) Ψ((Subscript[x, i]), (Subscript[y, i]), t)

The equation for the stream function is
	(∂_(x,x)) Ψ + (∂_(y,y))Ψ = (∑_(i=1)^N)(Subscript[Γ, i])δ(x - (Subscript[x, i]))δ(y - (Subscript[y, i])) ≈ ρ(x, y, t) 
where (Subscript[Γ, i]) is the vorticity of the ith line.

Numerics

For simplicity, we'll make the following assignments:
	Δx = Δy = Δt = 1 (the spatial and temporal grid points are integers)

With Δx = 1, the total length of our "box" of vorticies must be L >> 1.

The finite-difference approximation to the equations of motion are:
	Ψ(i+1,j) - 2Ψ(i,j) + Ψ(i-1,j) + Ψ(i,j+1) - 2Ψ(i,j) + Ψ(i,j-1) ≈ ρ(i,j)
	(Subscript[x, i])(t + 1) ≈ (Subscript[x, i])(t )  - (Ψ(i,j+1) - Ψ(i,j-1))/2
	(Subscript[y, i])(t + 1) ≈ (Subscript[y, i])(t )  + (Ψ(i+1,j) - Ψ(i-1,j))/2

In order to keep this solution numerically stable, the motion of the line vorticies cannot move too quickly. We’ll need to adjust the magnitude of Γ to keep the average motion sufficiently slow.

Basic Definitions

boxLength = 100; (* Length of plasma and number of grid points *)

numLines = 1000; (* equal numbers of positive and negative lines *)

lineΓ = 0.05; (* we need to adjust this to keep numerical errors small.. *)

SetDirectory[NotebookDirectory[]]

Negative Lines

Each particle is represented by a coordinate, {(Subscript[x, i]), (Subscript[v, i])}.

circleBlob[radius_,xCenter_,yCenter_] := Module[{r,ϕ},
r = RandomReal[{0,1}];
ϕ = RandomReal[{0,2π}];
{xCenter + radius r Cos[ϕ],yCenter + radius r Sin[ϕ]}]

negLines=Table[circleBlob[boxLength/10,boxLength/3, boxLength 2/3],{numLines}];
posLines=Table[circleBlob[boxLength/10,boxLength 2/3, boxLength 1/3],{numLines}];

Dimensions[negLines]

ListPlot[{posLines,negLines},AxesLabel->{"x","y"}, AspectRatio -> Automatic, Frame -> True, PlotStyle -> {Blue,Red}, PlotRange -> {{1,boxLength+1},{1,boxLength+1}}]

Vorticity Density

The vorticity density is evaluated on the grid points. We use "triangular weighting" to map the electron's position (in-between the grid points) to the charge density at each point.

rho0 = Table[0.0, {boxLength+1}, {boxLength + 1}]; (* including 0...boxLength *)
rhoT = rho0;

addρ[x_Real,y_Real] := Module[{locx, locy},
locx = Floor[x];
locy = Floor[y];
rhoT[[locx+1,locy+1]] += (locx+1-x)(locy+1-y);
rhoT[[locx+1,locy+2]] += (locx+1-x) (y - locy);
rhoT[[locx+2,locy+1]] += (x - locx)(locy+1-y);
rhoT[[locx+2,locy+2]] += (x - locx) (y - locy);]

getρ[vort_List] := Module[{},
rhoT = rho0;
addρ @@@ vort;
lineΓ rhoT]

(ρNeg = getρ[negLines]); // Timing

Plus @@ Flatten[ρNeg] (* How many lines? *)

ρPos = getρ[posLines];

ListContourPlot[ρPos-ρNeg, PlotLabel -> "Vorticity Density", ColorFunction -> "Rainbow", Contours -> 20]

(ρPos-ρNeg) // Dimensions

Linear Solver for Stream Function

The linear PDE for the stream function, Ψ(x, y), can be solved using the “Solve” command. Boundary conditions can be easily included when using a linear solve of the finite difference representation of the PDE. But, it is slightly slower, unless more creative algorthyms are used (e.g. multi-grid.)

eq[i_Integer,j_Integer] := (ψ[i+1,j] - 2 ψ[i,j] + ψ[i-1,j]) + (ψ[i,j+1] - 2 ψ[i,j] + ψ[i,j-1]) == ρ[i,j]

eq[3,4]

Boundary Conditions

ψ[1,_] = 0;
ψ[_,1] = 0;
ψ[boxLength+1,_] = 0;
ψ[_,boxLength+1] = 0;

Equations

eqns := Table[eq[i,j], {i,2,boxLength},{j,2,boxLength}]

unknowns = Flatten[Table[ψ[i,j],{i,2,boxLength},{j,2,boxLength}]];

Dimensions[unknowns]

unknowns // Short

Test Solution

rr = ρPos-ρNeg;

Do[ρ[i,j] = rr[[i,j]], {i,2,boxLength},{j,2,boxLength}];

es = Flatten[eqns];

Dimensions[es]

es[[{1,2,3,4}]]

sol = First[Solve[es,unknowns]]; // Timing

sol // Short

sol[[{110,111,112}]]

Plot

ψ[11,5] /. sol

ΦArray = Table[(ψ[i,j] /. sol), {i,1,boxLength+1},{j,1,boxLength+1}];

ΦArray[[11,5]]

Table[ΦArray[[i,20]], {i,1,10}]

ListPlot[Table[ΦArray[[i,20]], {i,1,boxLength+1}],PlotRange -> All]

ListContourPlot[ΦArray,PlotLabel->"Stream Function",AxesLabel->{"x","y"}]

FFT Solver

Another way to solve Poisson’s equation is to use the FFT. Since the FFT is a efficient algrorthym, this can speed the solution for Ψ, but we have to be clever to include boundary conditions. The easiest approach is “periodic” boundary conditions because the built-in FFT/IFFT algorythms can be used.
The inverse Fourier transform (Subscript[u, r]) of a list (Subscript[v, s]) of length n is defined to be (1/(Sqrt[n]))(∑_(s=1)^n)(Subscript[v, s])(e^(-2πi(r-1)(s-1)/n)).
The discrete Fourier transform (Subscript[v, s]) of a list (Subscript[u, r]) of length n is by default defined to be (1/(Sqrt[n]))(∑_(r=1)^n)(Subscript[u, r])(e^(2π i(r-1)(s-1)/n)).

kFourier2 = Table[-N[2 Sin[ π (nx-1)/(boxLength+1)]  + 2 Sin[π(ny-1)/(boxLength+1)]], {nx, 1, boxLength+1},{ny, 1, boxLength+1}]; (* This represents the "effective" wavenumber used to invert Poisson's equation with the FFT. *)

kFourier2 // Dimensions

Position[Flatten[kFourier2],i_/;i==0]

kFourier2[[1,1]] = 1.0;

Fourier[ρPos-ρNeg];//Timing

The FFT is ten times faster!

(ΦArray = Fourier[ρPos-ρNeg]/kFourier2;
ΦArray[[1,1]] = 0.0;
ΦArray = Re[InverseFourier[ΦArray]];)// Timing

ListPlot[Table[ΦArray[[i,20]], {i,1,boxLength+1}],PlotRange -> All]

ListContourPlot[ΦArray,PlotLabel->"Stream Function",AxesLabel->{"x","y"}]

Routine Parts

For a given vorticiyt density, this computes the electric field using an FFT to invert Poisson's equation...

getΨ[rhoPos_, rhoNeg_] := Module[{Φ},
Φ = Fourier[(rhoPos - rhoNeg)]/kFourier2;
Φ[[1,1]] = 0.0;
Φ = Re[InverseFourier[Φ]]]

(stream = getΨ[ρPos,ρNeg];) // Timing

Dimensions[stream]

getUx[stream_] := -(RotateLeft[stream,{0,1}] - RotateRight[stream,{0,1}])/2

getUy[stream_] := +(RotateLeft[stream,{1,0}] - RotateRight[stream,{1,0}])/2

sUx = getUx[stream];
sUy = getUy[stream];

vecs = Table[{sUx[[i,j]],sUy[[i,j]]}, {i, 1, boxLength+1}, {j,1,boxLength+1}];

ListVectorPlot[vecs, VectorPoints -> Fine, PlotLabel -> "Fluid Flow"]

Line Vortex Step

In this section, we use linear interpolation on a two-dimensional grid to find the fluid flow at any location.

stepLine[ux_List,uy_List][{x_,y_}] := Module[{xloc,xlocp, yloc,ylocp, vx, vy},
xloc = Floor[x];
xlocp = Mod[xloc+1,boxLength]+1;
yloc = Floor[y];
ylocp = Mod[yloc+1,boxLength]+1;
vx = ux[[xloc+1,yloc+1]] (xloc+1 - x)(yloc+1 - y) + ux[[xlocp,yloc+1]] (x - xloc)(yloc+1 - y)+ux[[xloc+1,ylocp]] (xloc+1 - x)(y- yloc) + ux[[xlocp,ylocp]] (x - xloc)(y- yloc) ; (* Interpolate for Ux *)
vy = uy[[xloc+1,yloc+1]] (xloc+1 - x)(yloc+1 - y) + uy[[xloc+1,ylocp]] (xloc+1 - x)(y- yloc) +  uy[[xlocp,yloc+1]] (x - xloc)(yloc+1 - y) + uy[[xlocp,ylocp]] (x - xloc)(y- yloc); (* Interpolate for Uy *)
(* Finally, advance position, periodically, … *)
{Mod[x + vx,boxLength,1], Mod[y + vy,boxLength,1]}];

negLines // Dimensions

(negLines1 = stepLine[sUx,sUy][#]& /@ negLines;) // Timing

(negLines1 = ParallelMap[stepLine[sUx,sUy][#]&, negLines];) // Timing

(posLines1 = stepLine[sUx,sUy][#]&  /@ posLines;) // Timing

negLines1-negLines // Dimensions

View of the change in the particles' phase-space positions...

ListPlot[{negLines1-negLines,posLines1-posLines} ,Frame->True, AspectRatio -> Automatic, PlotLabel -> "Single Time Step Motion"]

Many Steps

Now, we can repeat these calculations for many times steps. Be warned: as we increase the number of particles, this becomes a very long computation...

numSteps = 30;

pInd1 = 0;ProgressIndicator[Dynamic[pInd1],{1,numSteps}]
Timing[Do[
Do[ρPos=getρ[posLines];ρNeg=getρ[negLines];stream=getΨ[ρPos,ρNeg];
sUx = getUx[stream];
sUy = getUy[stream];negLines=stepLine[sUx,sUy][#]& /@ negLines;
posLines=stepLine[sUx,sUy][#]& /@ posLines;,{5}];
pInd1 = i;
lPlot[i] =  ListPlot[{posLines,negLines},AxesLabel->{"x","y"}, AspectRatio -> Automatic, Frame -> True, PlotStyle -> {Blue,Red}, PlotRange -> {{1,boxLength+1},{1,boxLength+1}}],{i, numSteps}];]

Animate[Show[lPlot[i]], {i, 1, numSteps, 1}, AnimationRunning -> False]

ListContourPlot[stream,PlotLabel->"Stream Function",AxesLabel->{"x","y"}]

ListContourPlot[ρPos-ρNeg, PlotLabel -> "Vorticity Density", ColorFunction -> "Rainbow", Contours -> 20]

Closer Interactions

What happens with a different initial condition?

In this section, we examine counter-drifting electrons.

halfcircleBlob[radius_,xCenter_,yCenter_,rl_:1] := Module[{r,ϕ},
r = RandomReal[{0,1}];
If[rl == 1,ϕ = RandomReal[{0,π}],ϕ = RandomReal[{0,-π}]];
{xCenter + radius r Cos[ϕ],yCenter + radius r Sin[ϕ]}]

negLines=Table[halfcircleBlob[boxLength/4,boxLength/2, boxLength /2,-1],{numLines}];
posLines=Table[halfcircleBlob[boxLength/4,boxLength /2, boxLength /2],{numLines}];

ListPlot[{posLines,negLines},AxesLabel->{"x","y"}, AspectRatio -> Automatic, Frame -> True, PlotStyle -> {Blue,Red}, PlotRange -> {{1,boxLength+1},{1,boxLength+1}}]

Many Steps

numSteps = 150;

pInd1 = 0;ProgressIndicator[Dynamic[pInd1],{1,numSteps}]
Timing[Do[
Do[ρPos=getρ[posLines];ρNeg=getρ[negLines];stream=getΨ[ρPos,ρNeg];
sUx = getUx[stream];
sUy = getUy[stream];negLines=stepLine[sUx,sUy][#]& /@ negLines;
posLines=stepLine[sUx,sUy][#]& /@ posLines;,{5}];
pInd1 = i;
lPlot[i] =  ListPlot[{posLines,negLines},AxesLabel->{"x","y"}, AspectRatio -> Automatic, Frame -> True, PlotStyle -> {Blue,Red}, PlotRange -> {{1,boxLength+1},{1,boxLength+1}}],{i, numSteps}];]

Animate[Show[lPlot[i]], {i, 1, numSteps, 1}, AnimationRunning -> False]

ListContourPlot[stream,PlotLabel->"Stream Function",AxesLabel->{"x","y"}]

ListContourPlot[ρPos-ρNeg, PlotLabel -> "Vorticity Density", ColorFunction -> "Rainbow", Contours -> 20]

Two Rings

What happens with a different initial condition?

In this section, we examine counter-drifting electrons.

ringBlob[radius_,xCenter_,yCenter_,thick_:2] := Module[{r,ϕ},
r = RandomReal[{-1,1}];
ϕ = RandomReal[{0,2π}];
{xCenter + (radius + thick r) Cos[ϕ],yCenter + (radius + thick r) Sin[ϕ]}]

negLines=Table[ringBlob[boxLength/3,boxLength/2, boxLength /2],{numLines}];
posLines=Table[ringBlob[boxLength/4,boxLength /2, boxLength /2],{numLines}];

ListPlot[{posLines,negLines},AxesLabel->{"x","y"}, AspectRatio -> Automatic, Frame -> True, PlotStyle -> {Blue,Red}, PlotRange -> {{1,boxLength+1},{1,boxLength+1}}]

Many Steps

numSteps = 150;

pInd1 = 0;ProgressIndicator[Dynamic[pInd1],{1,numSteps}]
Timing[Do[
Do[ρPos=getρ[posLines];ρNeg=getρ[negLines];stream=getΨ[ρPos,ρNeg];
sUx = getUx[stream];
sUy = getUy[stream];negLines=stepLine[sUx,sUy][#]& /@ negLines;
posLines=stepLine[sUx,sUy][#]& /@ posLines;,{5}];
pInd1 = i;
lPlot[i] =  ListPlot[{posLines,negLines},AxesLabel->{"x","y"}, AspectRatio -> Automatic, Frame -> True, PlotStyle -> {Blue,Red}, PlotRange -> {{1,boxLength+1},{1,boxLength+1}}],{i, numSteps}];]

Animate[Show[lPlot[i]], {i, 1, numSteps, 1}, AnimationRunning -> False]

(* lPlots = Table[Show[lPlot[i]], {i, 1, numSteps, 1}];
Export["ring_vortex.avi",lPlots]; *)

ListContourPlot[stream,PlotLabel->"Stream Function",AxesLabel->{"x","y"}]

ListContourPlot[ρPos-ρNeg, PlotLabel -> "Vorticity Density", ColorFunction -> "Rainbow", Contours -> 20]

Summary

“Despite subtle differences, the experiment and simula- tion show good quantitative agreement in several key areas.”